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SOLUTION SCHEME FOR TIME DEPENDENT 
HYDRODYNAMIC PLASMA FLOW ALONG 
A MAGNETIC FIELD LINE 
by 

J. M. Grebowsky 
R. K. Smith 

ABSTRACT 

A mathematical procedure, using the method of characteristics, is developed 
for solving the hydrodynamic flow equations in a nonhomogenous magnetic field 
for plasma flow along a field line in the presence of a gravitational field. Along 
a dipole magnetic field line with two supersonic plasma streams symmetrically 
directed towards the magnetic equator (a situation applicable to the earth's en- 
vironment) the solution is obtained by assuming the existence of shock disconti- 
nuities which propagate along the field line. The solution technique for the time 
development of the plasma flow is considered for both adiabatic and isothermal 
flow. A sample calculation for an adiabatic flow state shows explicitly the time 
evolution of adiabatic cooling and shock heating along a magnetic field line. 
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SOLUTION SCHEME FOR TIME DEPENDENT 
HYDRODYNAMIC PLASMA FLOW ALONG 
A MAGNETIC FIELD LINE 

Introduction 

Usually a sharp drop in the ambient electron density occurs with increasing 
altitude at a geocentric distance of a few earth radii in the earth's equatorial 
plane. This drop, referred to as the plasmapause, corresponds to the boundary 
separating plasma that drifts always across closed (i.e., both ends intersect the 
surface of the earth) magnetic field lines from plasma that at some time in its 
motion drifts onto polar magnetic field lines which connect the earth to the inter- 
planetary magnetic field (Nishida, 1966). On the latter field lines - the "open" 
lines - plasma is lost from the earth's environs by transport along the lines of 
force. This depletion along the open polar field lines occurs in the form of a 
supersonic polar wind of H^ ions directed away from the earth (Banks and Holzer, 
1968). 

Since the plasma which is transported along the open polar field lines also has a 
drift component perpendicular to the magnetic field lines due to the solar-wind 
magnetosphere interaction, the supersonic polar wind will convect onto closed 
magnetic field lines as shown in Figure 1. Assuming symmetry about the mag- 
netic equatorial plane, supersonic streams directed away from the earth in both 
hemispheres will be convected onto closed magnetic field lines. These streams 
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will collide near the equator producing collisionless plasma shocks (Banks, 
et al., 1970) which propagate down the field lines towards the earth. 

This paper will describe a mathematical scheme for solving the differential 
equations which govern the time development of plasma flow along a closed 
magnetic dipole field line when two supersonic streams directed along the field 
line are assumed initially to collide at the equator forming plasma shock waves. 
Although the differential equations to be considered are expressed in terms of 
a plasma state only, they also describe the flow of a compressible fluid along a 
channel with varying cross section. Hence the solution scheme to be described 
also has applicability within the field of fluid dynamics. 

System of Equations 

Only the plasma flow at high altitudes (i.e., above a few thousand kilometers) 
where the plasma consists of essentially ions and electrons is considered 
for simplicity. Chemical production and loss processes are ignorable in this 
region. If adiabatic flow is assumed (isothermal flow will be treated later) the 
hydromagnetic relations (see for example Montgomery and Tidman, 1964) 
describing plasma flow along the direction of the magnetic field can be written: 

+ vp^ + pv^ 

/?V^ + /3V V^ + P 3 = - ( 1 ) 
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where an isotropic pressure and a one dimensional coordinate system along the 
field line is assumed. The first of these equations is the continuity equation for 
flow in a tube (i.e., a magnetic flux tube) of varying cross section; the second 
is the momentum equation; and the last corresponds to the energy equation for 
adiabatic flow. The parameter t denotes time, s is the coordinate distance along 
the field line under consideration, p is the mass density (equal approximately to 
nM where n is the electron-ion number density and M is the proton mass), v is 
the flow velocity along the magnetic field direction, B is the magnetic field, P 
is the total plasma pressure (P = nkT where T is the sum of the electron and 
ion temperatures and k is Boltzman's constant), g is the component of the gravi- 
tational acceleration along the field line and y is the ratio of specific heats. The 
subscripts denote partial differentiation. 

Since the plasma parameters usually measured by satellite probes are in number 
densities and plasma temperatures, a more convenient form for these conser- 
vation relations is: 


n + vr»g + nVg = nvF 

kT k „ /ov 

nv^ + nvv^ + — +- nT^ = -ng (2) 

n(T^ + ) + (1 - y) T (n^ + vn^ ) = 0 

where F is defined as B^/B. This system of equations is a set of quasilinear 
partial differential equations of the first order for the three dependent variables 
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V, n and T which are functions of the two independent variables t and s. These 
equations are to be transformed into a system more suitable to a numerical 
solution. 

Characteristic Equations 

A system of hyperbolic differential equations may be transformed into a system 
of total differential equations by choosing appropriate paths in the s-t plane 
(see Courant and Friedrichs, 1948). As is shown in the Appendix, the hydro- 
dynamic equations (2) are hyperbolic provided, as the case in real physical sys- 
tems, the temperature remains positive. Hence the set of partial differential 
equations previously developed can be transformed into a system of total differ- 
ential equations along characteristic directions in the s-t plane. 

The characteristic directions for the hydrodynamic equations can be obtained 
directly by forming linear combinations of these equations and determining 
which of these combinations involve differentiation of the three dependent vari- 
ables V, n and T along the same direction (t (cr), s (cr)). For convenience the 
number density n will be expressed in terms of C = In n and each equation will 
be labeled as 


Lj = +vv^ 




k T 
M ' 


+ g = 0 


(3) 


Lj = - vF = 0 


(4) 
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^ = (S. + vc.)+ — +vT,) = o. 


(5) 


The linear combination can then be written 


L 







+ ^2) ''s + (^2 + ^3) + (^l T + ^2'' + 


(6) 


(i-r)T 


The equations describing the characteristic directions are, by inspection of (6) : 


Scr = 


^.2 + V 


X, / 


(7) 


X 


"" \Xj + X3 M 


1 kT _ 

TT- + V t. 


(8) 


Sa = 


^i(l-r) kT , 
: + t„. 

X, M ^ ^ 


( 9 ) 


Solving these equations for the ratios ^-2 Ai and X^/X^ when X j 0 yields the 
characteristic directions 


= V ± c 


(10) 
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where c s ykT/M is defined as the acoustic speed. The third characteristic 
direction corresponds to A. j = \ j = 0: 



( 11 ) 


For convenience these characteristic directions will be expressed as 


^ = V + /3c (12) 

^Or 

where ^ = +1, 0, -1. 

The corresponding linear combinations of the hydrodynamic equations (6) along 
the characteristic directions become: 


^ + :^ c i + (g - cvF) = 0 


(13) 


^ = 0 : 


(7-l)c " 


c_ = 0 


(14) 


/3=-l: i c 1 + (g + cvF) = 0 (16) 

where the subscript cr denotes differentiation along the characteristic,. The six 
characteristic equations - (12) through (15) - replace the original three flow 
equations (3) through (5). 
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Solution Scheme and Boundary Conditions 


Assume that the variables n, v and T are specified along a given boundary 
curve (this curve cannot be a characteristic) in s-t space. Consider a point 
(s, t) near the boundary curve (see Figure 2). If no discontinuities exist at 
this point, it will correspond to the intersection of three (i.e., /3 = -1, 0, 1) dis- 
tinct characteristics which also intersect the boundary curve. The variation of 
the variables n, v and T along these characteristic directions will be described 
by equations (13), (14) and (15). 

If the computation point (s, t) is selected very near the boundary curve, then the 
characteristics will be essentially straight lines between the boundary curve 
and (s, t) and will be described by ds/dt = v^ + /3 c^ where /3 = + 1, 0, -1. The 
subscripts now denote boundary values and not differentiation (e.g., Vj is the 
velocity at the point where the boundary curve is intersected by the /3 = 1 
characteristic passing through (s, t). Differentiation will be denoted explicitly 
in the work to follow so this change of notation should cause no difficulties. 

For each selected point (s, t) near the boundary, the boundary curve is searched 
to locate the characteristics which pass tlirough (s, t) and to determine the 
boundary values (i.e., n^, v^, T^ where (i = -1, 0, 1). The finite difference 
equations corresponding to the characteristic equations (13), (14) and (15) can then 
be solved at the point (s, t). The finite difference equations can best be solved 
if considered in matrix form: 


7 


1 


V 


b. 

y 7 




1 

0 1 _ ^ 




t>o 

1 

o 

o 




2 






j C_i o 


c 


b- 

7 7 
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where 

= vj + i Cl + - (ii - CjVjFj) Atj 

K r 2 

- ^0 - (y:rr) 

^3 = ''-1 - ■ 7 ‘'-1 ■ ‘'-1 ''-1 ^ 

and g^, F^ denote averages along the segment of the /i characteristic between 
the point (s, t) and the boundary whereas ^t^is the time increment along this 
characteristic segment. 

After the solutions are generated along a curve located a short distance from the 
initial boundary (e.g., along the dashed curve in Figure 2), then this curve can 
be treated as a new boundary curve for further calculations. By proceeding in 
this manner, the solution can be generated at each point in s-t space at which 
the three intersecting characteristics can be traced back to the initial boundary 
curve. It remains now to determine what type of boundary curve is best suitable 
for the problem under consideration. 
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Since the plasma state along the earth’s outer magnetic field lines will depend 
upon the rate at which ionization is supplied from the lower ionosphere and upon the 
rate at which the ionospheric parameters change, a natural choice for a boundary 
condition is the specification of n, v and T as a function of time at a base level 
altitude. This base level point (s^^) must be within the protonosphere since the 
flow equations have been developed assuming a neutral proton-electron plasma. 
However if the plasma flow is everywhere supersonic along the field line, this 
boundary condition alone is not sufficient to determine the solution at every 
point along the field line. This is due to the fact that the only points (at time t) 
which correspond to the intersection of three distinct characteristics emanating 
from the boundary curve s = s^ are those points with s coordinates smaller 
than or equal to the s position (at time t) of the yS = -1 characteristic which 
passes through the point (s^ , t = 0). As seen in Figure 3, these points (region 
C) do not cover all of s - t space. 

If, on the other hand, boundary conditions are also specified along the boundary 
t = 0, then the solution can be obtained everywhere in s - t space. This is 
readily seen in Figure 3, where region A depicts those points at which the solu- 
tion is obtainable using characteristics which cross the t = 0 boundary only, and 
where B defines the s - t region in which solutions are generated using a mixture 
of characteristics, some of which cross the t = 0 boundary and some of which 
cross the s = s^ boundary curve. 
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It should be noted that the t = 0 and s = s q boundary values must be selected 
so that they are compatible with one another. For example if subsonic flow 
exists near the point (s^, t = 0), then /3 = - 1 characteristics (i.e., ds/dt - v - c) 
emanating from the boundary curve segment t = 0 near (s^, o) may cross the 
curve s = s^. In such a case the value of n, v and T cannot be selected arbi- 
trarily on both boundary curve segments, but must be chosen so that equation 
(20) is satisfied. Even in the supersonic regime, if the boundary values are not 
selected properly, accidental discontinuities can be introduced near (s^, 0). 

From the physical viewpoint the most convenient computation scheme consists 
of determining the plasma state at equally spaced points along the entire field 
line as a function of time. That is, initially the solution is obtained at all of the 
selected points along the field line at time At, where At is a time increment 
selected small enough that the characteristic segments between the initial 
boimdary t = 0 and the computation curve t = At in (s, t) space are approximately 
straight lines. The values of v, n and T determined at t = At are then used to 
generate the solution along the field line at t = 2 At, etc. This procedure will 
develop explicitly the time evolution of the plasma flow along the field line. 

Shock Problem 

As indicated in the introduction, the two supersonic plasma streams (one in the 
northern hemisphere, the other in the southern hemisphere) directed along the 
chosen field line away from the earth will interact initially at the equator 
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forming shock discontinuities. Assuming equinoctal conditions (i.e., reflection 
symmetry of the plasma state about the magnetic equatorial plane) two shocks, 
one on each side of the equator, will be formed initially at infinitesimal distances 
from the equatorial plane. The existence of these shocks satisfies the symmetry 
condition at the equator (i.e., v = 0) as the upward directed supersonic flow be- 
comes subsonic in its passage through the shock. After the shocks have been 
formed initially at the equator, they will propagate along the field line towards 
the earth as shown in Figure 4. In the following sections only one of the two 
shocks will be treated explicitly - i.e., the plasma flow will be developed ex- 
plicitly along only half of the considered field line. The flow with its associated 
shock on the other half of the field line can be obtained by a simple mirror re- 
flection about the equatorial plane. 

Shock Relations 

In order to relate the plasma state on the earth side (i.e., upstream) of the shock 
to the plasma state on the equator side (i.e., downstream) use is to be made of 
the fact that conservation of the particle number flux, momentum flux, and 
energy flux must be maintained through the shock. If the plasma parameters 
on the earth side of the shock discontinuity are denoted by the subscript 0 and 
the equator side values by the subscript 1, these conservation relations are 
(see Courant and Friedricks, 1948, for the deviation of these relations from 
the fluid equations): 
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Pl^[ =Po^O 

(^iVi)v; - (PoVo)v' =Pq -Pj (17) 



where the primes denote flow velocity with respect to the shock which is moving 
with velocity (e.g., v ' = Vj - V^) and where e is the internal energy/mass of 
the plasma. It should be noted that the velocities are considered to be positive 
when directed along the field line away from the earth, 

hi equations (17) if 7 1 (the isothermal case 7=1 will be discussed in a later 

section) the energy density e can be expressed as 


1 kT ^ 1 

(7-1) M “ 7(7 - 1) ^ 


(18) 


whereas the pressure P, as defined earlier, is just nkT and the mass density P 
is nM. Using these expressions, equations (17) are expressible in terms of n 
and c rather than p , e and P: 




OiVjV; -n^VoV' = 1 (noc2 -n^c^) 




1 


2 ^ 7(7-1) * 


c? 1 - + 


(19) 


2 0 7(7 - 1) 


2 _ 
Co - 


— (noVocJ -niVjcJ). 
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When upward directed supersonic flow exists everywhere downstream of the 
shock, the solution to the flow equations on the earth side can be determined 
completely at all points using the previously discussed procedure. The shock 
discontinuity does not enter into these computations except as an upper cutoff 
boundary. However the plasma state on the equator side of the shock cannot be 
determined without considering the effects of the shock since all of the charac- 
teristics connecting this region to the s = s^ boundary would pass through the 
shock. In fact, as will be seen, the characteristics on the equator side of the 
shock are determined in part by using the shock path in s - t space as a boimdary 
curve. The corresponding boundary values on the equator side of the shock can 
be determined by solving the shock jump conditions - equations (19). 


Since the earth side flow can be determined without taking into account the 
shock whereas the equator side flow cannot, equations (17) are best solved for 
the parameters on the equator side of the shock in terms of the parameters up- 
stream of the shock. Such a solution is given by the following expressions: 


_"_o _ (y-i)vo + 2cg 

"i (y + 1) v'^ 


( 20 ) 


v' "i 


■I =(“) -'y^o - v^)] 


(21) 


(22) 


These equations correspond to the Rankine Hugoniot equations. 
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Equations (20) through (22) determine the downstream parameters at the shock 
only when the shock velocity or one other downstream parameter is known. Since 
a nonsteady state is under consideration, these equations will not, in general, 
form a closed system except at t = 0 when the shock is only an infinitesimal 
distance from the equator and v i = 0 by symmetry. In such a case equation (21) 
can be replaced by the relation for the shock velocity which is: 

= -a - /x2 +/32 (23) 

where 

a = j(7- 3 )Vq 
and 

/3= ^ (r- 1) + 2 c2 . 

When t 0 the downstream flow equations must be utilized to determine the con- 
ditions at the shock. This procedure will be discussed in the next section. 

Shock Solution 

Since the flow on the equator side of the shock depends upon the shock conditions 
and since the shock jump conditions do not in general form a closed set of equa- 
tions, another equation must be set up, downstream of the shock, which does not 
increase the number of dependent variables . One of the previously developed 
characteristic relations - equation (15) - describing the flow downstream of 
the shock will be used for this purpose. 
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Consider the situation near the shock as depicted in Figure 5. The flow at all 
points on the earth side of the shock, as demonstrated previously, can be de- 
termined without reference to the shock. Now assume that the solution has been 
obtained at all points on the field line downstream of the shock at time t Then 
since the flow is subsonic at these locations, /3 = - 1 characteristics connect 
downstream points at t = t^ to the shock path between t^ and t^ + A t (see Figure 
5). Hence equation (15) des ribing variations along these characteristics can be 
used to complete the solution to the jump conditions. The /3 = + 1, 0 character- 
istics, on the other hand, have positive slopes in s - t space and do not connect 
the t = t^ downstream points to the shock and hence do not directly relate to the 
solution for the shock motion. 

Assume that the shock position at time t^, denoted by X^, and the corresponding 
shock velocity are known (i.e., they were previously computed). Then if the 
shock velocity changes only slightly in the time interval A t, the shock velocity 
at t = tp + At can, to a good approximation, be expressed as a linear function of 
which denotes the position (as yet unknown) of the shock at time t^ + At. 

That is: 

2(Xj -Xq) 

Vs(Xi) - It • -^s(Xo) (24) 
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Since the flow solution is Icnown everywhere on the earth side of the shock, the 
variables Uq, Vq, Cq are determinable functions of s and hence of Xi at time 
tp + A t. Substituting equation (24) into the shock relations (20) through (22) will 
then yield three simultaneous equations for the four unknowns nj, Vj, c^ and X j. 


The dependent variables along the field line on the equator side of the shock at 
t = tj. are assumed to have been previously determined and are thus laiovm 
functions of the distance along the field line - in this region these variables are 
denoted by the subscript I. Therefore the /^ = - 1 characteristic passing through 
the shock point (X^, t + A t) relates the downstream shock parameters n j, v^, c^ 
and the shock position X^ to the position s ^ at which this characteristic crosses 
the t = t line: 

c 


S 


1 




(25) 


where Vt, Cj are functions of Sj, The finite difference form of equation (15) 
along this characteristic segment is; 


V,- Vj . _ 


1 *'1'^ I 2 I 

(n,- Oj) -_(c - c ) + lg + 

7 Oj+ n j * * 7 * ^ ' 


(c,+ c,) (v,+ v,) 

F I At = 0. (26) 


This equation in conjunction with equation (25) yields a relation between Vj , n^ 
c^ and the shock position X Hence these relations and the shock jump condi- 
tions form a closed set of equations which can be solved for the shock position 
and the downstream shock variables. 
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Downstream Solution 


The solution at points between the shock and the equator can readily be obtained 
using the method of characteristics described earlier. However, the boundary 
curve is not the same as was used for the upstream solution. The t = t^ line 
segment in s *■ t space on the equator side of the shock is still a valid boundary 
curve segment for computing the solution at a later time t = t,. + ^ t. However 
the s = s^ boundary curve segment used previously is not directly applicable to 
the downstream solution. It is replaced by the shock path in s - t space. Hence 
the shock position at time t^ + A t must be computed along with the equatorward 
shock variables before the solution can be obtained downstream from the shock. 
The variation of n^, v^, c^ along the shock path is obtained by interpolation be- 
tween the values at t„ and t„ + At. 

C C 

Isothe rmal Flow 

The flow solution thus far has been developed explicitly for the case 7 ^ 1 * If 
an isothermal state (y = 1) exists along the field line (ignoring the shock mo- 
mentarily) the solution scheme must be modified slightly. When y = 1 the 
temperature, which is specified as a function of time at the base level s^(i.e., 
along the s = Sq boundary curve segment), is known at all points along the field 
line and only the /3 = + 1 and -1 characteristic equations are required for a 
solution. The energy equation (14) is no longer applicable nor is it needed for 
a solution. Otherwise the solution is obtained in the manner previously 
described. 
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when the shock discontinuities are included and isothermal conditions prevail 
upstream and downstream of the shocks, the field line is separated into two dis- 
tinct isothermal regions each with its own characteristic temperature. The 
spatially invariant temperature on the equator side of the shock is coupled to the 
upstream temperature via the shock jump conditions. The shock process, how- 
ever, is not isothermal and hence cannot be described by simply setting 7 - 1 
in the shock relations (19). The conservation conditions (17) are still valid 
when the temperature is spatially invariant on both sides of the shock, and the 
first part of (18) relating the internal energy to the temperature is correct if 
y is taken to be the value leading to the internal energy per mass (e) of the 
plasma (e.g,, if three internal degrees of freedom exist y = 6/3), However the 
relationship between e and the sound speed c must be modified since the iso- 
thermal sound speed on either side of the shock is c ^ = kT/M and hence is not 
related to the factor y describing the energy transfer through the shock in the 
first equality of equation (18), That is, (18) is replaced by 

e = ^ 

(r- 1 ) 

where 7 / 1. When this change is made, the modified form of the jump conditions 
( 20 ) through ( 22 ) are readily obtainable and the solution scheme then follows the 
procedure described previously using the /3 = +1 and -1 characteristics. 
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Sample Computation 


A sample calculation will show explicitly the type of time evolution which is 
sought. The solution is developed along a field line of the earth’s magnetic 
field (assumed dipolar) which intersects the equatorial plane at a geocentric 
distance of 6 earth radii (1 earth radius corresponds to 6370 kilometers). For 
simplicity the dependent variables n, v and T are fixed at constant values along 
the boundary curves (i.e., along the field line initially, and at all times at the 
base level- selected to be at an altitude of 3000 km). The fixed values selected 
were: 

n = 300 particles /cc 

V - 11 km/sec 

T = 8000°K 

These values were selected arbitrarily and are not to be taken as representative 
of any real geophysical state. 

With these boundary conditions the flow develops along the field line as shown 
in Figure 6. The shock which is formed initially at the equator is accelerated 
towards the base of the field line leaving behind a hot, relatively dense, subsonic 
flow. In this case a steepening compression wave is seen to propagate towards 
the shock. The computed density and temperature on the earth side of the shock 
in a short period of time approach a state of rapid decrease with increasing 
distance along the field line. This results from the adiabatic expansion of the 
plasma along the diverging magnetic flux tube. These computations indicate the 
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type of flow behavior which can be explored using the solution scheme discussed 


in this paper. 


APPENDIX A 


The equations under coneide ration can be expressed in the general form 

L = u<^^> + f 0 (A-1) 

where the summation convention is assumed ( cr is summed over the number of 
dependent variables, v is summed over the number of independent variables). 
The greek letters denote labeling indices whereas the subscript denotes 
differentiation with respect to the independent coordinate • The coefficients 
%(T general, functions of the independent variable and the dependent 

variable u 

In equations (A-1) the functions u^°"^ are differentiated in different directions. 

It would simplifj' the numerical computations if characteristic surfaces could 
be determined for which a linear combination L = of these equations in- 

volves derivatives of the functions u only in the direction of the surface ele- 
ment. The condition for the existence of the multipliers is the determinant 
equation 

II II =0 (A-2) 

for the components of the normal vector { ^ characterizing the surface ele- 
ments (see Courant and Friedricks, 1948). 


The system of equations (A-1) is identical to equations (3), (4) and (5) if the vari- 
ables are defined as: 


/Lt = 1, 2, 3; (7 = 1, 2, 3; u = 1, 2 


( = v(s, t) 

f , = ng 

^2) = n(s, f) 

fj = -nvF 

=T(s, t) 

II 

O 


(A-3) 


and the matrix elements a^^ where denotes the row and a the column are 
given by 


and 


a 


( 1 ) 

tier 


0 1 0 

u<2) 0 0 

P (1 - 


(A-4) 


|( 2 ) 


u 


(i) 


a(2) = 
fia 


u(2)|j(l) ii u (3) 
M 


^ U<2) 

M 


0 


(1 - y)u ^ ^ ^ u^2)jj(i) 


Substituting these relations into equation (A- 2) and expanding the determinant 
yields the relation 


( u ( 2))2 ^ ) 


7 — _ 7^2 _ 7 ] + (U^^^)2 

M 


= 0 


(A-5) 
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where 


4 ( 2 ) 

defines the direction of the vector { ^ . 


One solution of (A-5) is readily seen to be t] = -u which is always real. The 
other solutions are obtained from the quadratic relation formed when the square 
bracketed term in (A-5) is set to zero. Both solutions of the quadratic equation 
in 7 ] are real if yk/M u 0. Since u<3) = T and since the temperature is 
always positive in a real physical state, this latter relation is always satisfied. 
Hence three distinct characteristic surface elements exist and the system of 
hydrodynamic equations imder consideration are totally hyperbolic (Courant and 
Hilbert, 1937), making a numerical solution along characteristic directions possible. 
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Figure 3. Region A Is Covered by Characteristics Emerging from the t a to 
Boundary, Region B Is Connected to both the t = t© and s = Sq 
B oundary Curve Segments, Whereas Region C Is Coupled to 
Only the s = Sq Boundary. The Subscript b Denotes Values 
Along the Characteristic Emerging from (so,to). 
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Figure 4. The Shocks Produced Initially at the Equator Propagate 
Towards the Earth Leaving Behind a Subsonic Flow, 
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Figure 5. The Figure Shows That /J = -1 Characteristics (clashed lines) 
Connect Computed Points (x's) to the Shock Path Whereas 
the |3 = +1 Characteristics (dotted curves) on the Equator 
Side of the Shock do not. 
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Figure 6. Flow Development Along a Dipole Field Line Which Crosses the Earth's Equatorial Plane 
at 6 Earth Radii with Constant Dependent Variables Along the Boundary Curves 
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